function heat_capacity_versus_T()
    n = 40;
    p0 = 0.0;
    b = 0;
    isingmat0 = isingmat_ini(n,p0);
    flip_num = 1;
    %img_num = 3;
    
    steps = 20000;
    sam_len = 8000;
    col_step = 5000;
    
    num = 11;
    T_min = 2.3;
    T_max = 2.5;
    
    T_array = linspace(T_min, T_max, num);
    
    Clist = zeros(1,num);
    
    mat_list = cell(1,num);
    
    for i = 1:num
        isingmat = isingmat0;
        T = T_array(i);
        for cy = 1:steps
            [trans_mat, km, mm]= trans_state(isingmat, n, T, b, flip_num);
            isingmat = trans_mat;
        end
        
        mat_list{i} = isingmat;
        
        en_list = zeros(1,sam_len);
        
        isingmat_T = isingmat;
        
        if(T<1.8)
            re_col_step = 300;
        else
            re_col_step = col_step;
        end
        
        for p = 1:sam_len
            isingmat = isingmat_T;
            for q = 1:re_col_step
                [trans_mat, km, mm]= trans_state(isingmat, n, T, b, flip_num);
                isingmat = trans_mat;
            end
            en = ising_energy(isingmat, n, b);
            en_list(p) = en;
        end
        
        C = var(en_list)/(T*T); 
        Clist(i) = C;
    end
    
    %将mat_list写成txt文件
    fname = sprintf('p_%.1f-inisteps_%d-samlen_%d-colstep_%d.txt', p0, steps, sam_len, col_step);
    fid = fopen(['.\',fname],'w');
    for k = 1:num
        mat = mat_list{k};
        for ii = 1:n
            for jj = 1:n
                fprintf(fid, '%d\t',mat(ii,jj));
            end
            fprintf(fid, '\n');
        end
    end
    fclose(fid);
    
    %将T_array和Clist写成txt文件
    fname1 = sprintf('C(T)_data--p_%.1f-inisteps_%d-samlen_%d-colstep_%d.txt', p0, steps, sam_len, col_step);
    fid1 = fopen(['.\',fname1],'w');
    for i = 1:num
        fprintf(fid1, '%d\t',T_array(i));
    end
    fprintf(fid1, '\n');
    for i = 1:num
        fprintf(fid1, '%d\t',Clist(i));
    end
    fprintf(fid1, '\n');
    fclose(fid1);
    
    figure(1)
    scatter(T_array, Clist)
    xlabel('T');
    ylabel('C (a.u.)');
    title_string = sprintf ('C(T) curve (n = %d, b = %.1f, p = %.1f)', n, b, p0);
    title(title_string);
    
%     img_step = round((num-3)/(img_num-1));
%     for i = 2 : img_step : (num-1)   
%         figure(i)
%         ising_plot(n, mat_list{i});
%         title_string = sprintf ('configuration of T %.1f (n = %d, b = %.1f)', T_array(i), n, b);
%         title(title_string);
%         pause(0.5);
%     end

    end

    
    
    